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ABSTRACT 



The equation of state (EOS) of the solar interior is accurately and smoothly 
determined from ab initio simulations named quantum Langevin molecular dy- 
namics (QLMD) in the pressure range of 58 < P < 4.6 x 10 5 Mbar at the 
temperature range of 1 < T < 1500 eV. The central pressure is calculated, and 
compared with other models. The effect of heavy elements such as carbon and 
oxygen on the EOS is also discussed. 

Subject headings: equation of state — solar interior — ab initio 
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1. Introduction 

For solar and stellar models, a high-quality equation of state (EOS) is very crucial 
(Basu & Antia 2008). It is well-known that some requirements should be satisfied for solar 
and stellar modeling: thermodynamic quantities would be smooth, consistent, valid over a 
large range of temperature and density, and incorporate the most important astrophysically 
relevant chemical elements (Dappen 2006) as well. There are two major efforts for a 
high-precision and high-accuracy EOS made before, which have been included in the recent 
opacity recalculations. They are the international Opacity Project (OP) (Seatom 1995; 
Berrington 1997) and Opacity Project at Livermore (OPAL). In OP, the model named 
Mihalas-Hummer-Dappen (MHD) (Hummer & Mihalas 1988; Mihalas et al. 1988; Dappen 
et al. 1988; Nayfonov et al. 1999; Trampedach et al. 2006) for EOS is developed, dealing 
with heuristic concepts about the modification of atoms and ions in a plasma. In the other 
one, OPAL, the EOS relies on physical picture, which is built on the modeling of electrons 
and nuclei. This model is called ACTEX (activity expansion) EOS (Rogers 1986; Rogers 
et al. 1996; Rogers & Nayfonov 2002; Iglesias & Rogers 1991). Both MHD and OPAL are 
dependent on the potentials between particles (electron-electron, ion-electron, ion-ion). In 
addition, there are other models such as Eggleton, Faulkner, & Flannery (1973) (EFF), 
which is thermodynamically consistent and qualitatively correct. EFF does not consider 
excited states, H molecule formation, or Coulomb corrections, and treats full ionization 
for heavy elements at high density. It is interesting to note that the effect of partially 
degenerate electrons can be included partly according to the Fermi-Dirac statistics. In most 
cases, the EOS from OP, OPAL and even EFF should be good for the input parameters 
of stellar models. However, there is always some physics such as the coupling of ions 
which is missed in these models. With the increasing requirements of the high-precision in 
helioseismic study, a parameter-free model beyond Debye-Hiickel approximation for a more 
accurate EOS than that of OPAL and MHD is still necessary (Dappen 2006; Dappen & 
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Nayfonov 2000). 

The conditions in the solar interior are very complicated, where the densities are from 
10 g/cm 3 up to 160 g/cm 3 , and temperatures from 50 eV to 1400 eV (Bahcall et al. 2001). 
At the same time, Hydrogen and Helium are the main elements, taking about 98%, with 
small abundances of other heavy elements such as carbon, oxygen and iron. In order to 
obtain an accurate and smooth EOS for the whole sun, it is necessary to develop a model 
which can cover all conditions in the sun. It is very clear that matter in the sun is not 
always ideal ionized gas plasma, but moderately coupled, partly degenerate, and partly 
ionized (especially for heavy elements) in some area according to the definition of coupling 
parameter T and degenerate parameter 9 (Ichimaru 1982), where V = Z* 2 /(ksTa), with 
T the system temperature, kg the Boltzmann constant, a the mean ionic sphere radius 
defined as a = (3/(47m i )) 1 / 3 , Z* the average ionization degree, rij the ionic number density, 
and 9 = T/Tp, with the Fermi temperature Tp = (37r 2 n e ) 2 / 3 /2 (n e is the number density 
of electrons). Generally speaking, when the degenerate parameter 9 ~ 1 and coupled 
parameter T ~ 1, the matter is considered to be partially degenerate and moderately 
coupled. Otherwise, if 9 <C 1 (9 ^> 1) and T ^> 1 {V 1), matter is strongly (weakly) 
degenerate and strongly (weakly) coupled. For most of regimes in solar interior, matter 
is weakly coupled and weakly degenerate. However, the extremely high density can not 
promise the negligibility of the coupling of ions, and the theory of ideal ionized gas plasma 
is not always appropriate. For weakly coupled or weakly degenerate matter, single atomic 
models such as average atoms (AA) models (Hou et al. 2006; Yuan 2002), and detailed level 
accounting (DLA) (Zeng & Yuan 2004) are assumed valid. Generally, AA model is built for 
the dense plasma, and DLA model is available for the relatively low density. Very recently, 
AA model was applied for the astrophysical plasmas successfully (Faussurier et al. 2010). 
However, there is no direct physical evidence supporting these approaches, and we don't 
even know whether they are correct in such dense matter, and how they behave for the 



-5 - 



extreme dense matter under conditions such as the solar center. Therefore, it is extremely 
important to develop a more reliable model for all these conditions. 

For the strongly (moderately) coupled matter with partly degenerate electrons, 
quantum molecular dynamics (QMD), without requiring any assumptions about the 
potential between atoms, supplies a powerful and accurate tool, and has been successfully 
applied in warm dense matter (WDM) (Collins et al. 1995; Desjarlais et al. 2002; Mazevet 
et al. 2005, 2008). In astrophysics, QMD has been used to study the properties of giant 
planets and given rise to amazingly satisfying results (Lorenzen et al. 2009; Militzer & 
Hubbard 2009; Gillan et al. 2006; Vorberger et al. 2007; Nettelmann et al. 2008; Wilson & 
Militzer 2010). Very recently, QMD was extended to the field of high energy density physics 
(HEDP) by considering the electron-ion collision induced friction (EI-CIF) in the procedure 
of molecular dynamics, and the corresponding model called quantum Langevin molecular 
dynamics (QLMD) was constructed (Dai & Yuan 2009; Dai et al. 2010). Thanks to this 
model, one can go into the regime of HEDP from first principles, and the thermodynamic 
properties of the sun and sun-like stars can be therefore researched based on ab initio 
method. It also gives us a good tool to investigate the system in which the coupling of ions 
is important. 

In this work, we firstly study the EOS of conditions under the solar interior using 
QLMD. Three densities with different temperatures are chosen, and their EOS are compared 
with the EOS of the AA model with energy-level broadening (AAB) (Hou et al. 2006). 
QLMD can ensure much more accurate results in all conditions, and agrees with AAB at 
high temperatures. The central pressure of the sun is also calculated using different models, 
comparing with the data of other models and standard solar models (Bahcall et al. 2001). 
The important effect of heavy elements on EOS of very dense matter is also investigated at 
the end. 
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2. 



Theoretical method 



2.1. Quantum Langevin molecular dynamics 



We firstly recall the QLMD method here, where the electronic structure is studied 
from density functional theory (DFT), and the ionic trajectory is performed using Langevin 
Equation (Dai et al. 2010) 



Where Mj is the ionic mass, F is the force calculated in DFT, 7 is a Langevin friction 
coefficient, Ri is the position of ions, and Ni is a Gaussian random noise corresponding to 
7. Considering the dynamical collisions between electrons and ions at high temperature, 
the friction coefficient 7 can be estimated by 



Where m e is the electronic mass, Z*, rii and T are the same as the above definition. 
Therefore, the dynamical electron-ion collisions at high temperature can be described 
as a friction or noise effect. By introducing this dynamical EI-CIF, the first principles 
simulations can be applied in the HEDP regime, overcoming the difficulty of numerical 
calculations, and therefore give more reliable results. 

To guarantee an accurate sampling of the Maxwell-Boltzmann distribution, the noise 
has to obey the fluctuation-dissipation theorem: 



Where dt is the time step of molecular dynamics. At the same time, the random forces are 
taken from a Gaussian distribution of mean zero and variance of < Ni 2 >= ^M^BTjdt. 

In order to integrate Eq. 1, the formalism in a Verlet-like form (Pastor et al. 1988) 



M/R,! = F - 7M 7 R! + Ni, 



(1) 




(2) 



< Ni(0)Ni(t) >= 6jM T k B Tdt 



(3) 
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integration is performed as follows 




(4) 



+ (dt 2 /Mj)(F BO (t) + N x (t)(l + )p T dt)-\ 



The velocities of ions vi(t + dt) can also be calculated by using the Verlet formula 



This ab initio molecular dynamics model based on LE is named QLMD, which extends 
the applications of ab initio method into the field of HEDP. Based on QLMD, the effects 
of coupled ions and degenerate electrons can be studied, and give much more accurate 
results at relatively low temperature and high density. It is worth pointing out that the 
computational cost of the QLMD is very expensive under conditions of very weakly coupled 
and non-degenerate matter. In this case, the existing models such as MHD, ACTEX, and 
AA model can give consistent and accurate data. 

2.2. Average atoms model with electronic energy-level broadening 

The average atom model is one of the statistical approximations applied to study the 
electronic structure of atoms and ions in hot dense plasmas, which is easily to be applied in 
conjunction with a variety of treatments of electron orbitals in atoms. In a full relativistic 
self-consistent field-based AA model, the influence of the environment on the atom is 
assumed to be spherically symmetric on average. The movement of an electron under the 
interactions of the nucleus and other electrons is approximated by a central field, which is 
determined by the standard self-consistent calculation. In the central field, the radial part 



Vl (t + dt) = Ri 



Ri(t + dt) - Ri(t - dt) 



(5) 



2dt 
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of the Dirac equation has the form: 

^ + ^(r) = l[e + c 2 - V(r)}Q nK (r) 

(6) 

^ - *Qn*(r) = -l[e - c 2 - V(r)]P nK (r) 
where P(r) and Q(r) are respectively the large and small components of the wave function, 
c is the light speed, and V(r) is the self-consistent potential, consisting of three parts, which 
are respectively the static, exchange and correlation potentials. The static part is calculated 
from the charge distributions in the atom, while the exchange and correlation parts take 
the approximate forms of Dharma-Wardana and Taylor (Dharma-Wardana & Taylor 1981). 
For bound states, we have the boundary conditions satisfied by the radial wave functions 

P nK (r) ^ ar* 1 f P nK {r) ^ ar^ 

or < (7) 

where Rj, is the radius of the atom. The electron distribution is calculated separately for 
the bound and free electron parts. The bound electron density is obtained according to 



h = — ' ^ ■ 7 - (9) 



b 

where bj is the occupation number of the state j. In A A model without energy- level 
broadening, the occupation number bj is determined by the Fermi-Dirac distribution 

exp((ej -fi)/T) + l' 

In order to consider the effect of energy-level broadening, Gaussian functions p(e) centered 
at the corresponding electron orbital energies of Eq. 7 are introduced into the Fermi-Dirac 
distribution of electrons, i.e., 

W = „ 2|Kjl ^\ (io) 



and 



= i^Ejf b M(P?(r) + Q^de, (11) 
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where 1 = J p(e)de and the two orbital energies obtained from the boundary conditions 
of Eq. 7 are taken as the upper and lower half maximum positions of the Gaussian form 
energy-level broadening. 

Based on this approach, the splitting of the real energy levels approximated by the 
energy-level broadening can be considered, which makes the irregularities caused by the 
pressure-induced electron ionization without energy-level broadening disappear naturally 
(Hou et al. 2006, 2007). 



3. Results and discussions 

3.1. Computational details 

For the study of solar interior EOS, we choose three typical densities in the solar 
interior: 10 g/cm 3 (radiative zone), 100 g/cm 3 and 160 g/cm 3 (core). Temperatures from 
1 eV to 550 eV for 10 g/cm 3 , and 10 eV to 1500 eV for 100 and 160 g/cm 3 are calculated, 
covering the conditions from the radiative zone into the core. Since the H and He elements 
take up about 98% of the composition of the sun (Lodders 2003), we calculated four 
structures with different compositions (using X to represent the mass abundance of H, 
and Y for abundance of He, Z for others), which are X=l, Y=0, Z=0; X=0, Y=l, Z=0; 
X=0.7, Y=0.3, Z=0; X=0.40, Y=0.60, Z=0. The supercells containing T=125 particles for 
10 g/cm 3 and T=256 particles for 160 g/cm 3 are constructed. When Z=0, the number of H 
atoms (Nh) and He (Nne) can De calculated by the formulas: 

-£f- = Y , NHe + NH = T (12) 

In these conditions, the matter goes from strongly coupled to relatively weakly coupled, 
and from partially degenerate to weakly degenerate, as shown in Fig. 1, where QLMD has 
been proved successfully. 
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For the QLMD simulations, we have used the Quantum Espresso package (Giannozzi et 
al. 2009) based on the finite temperature DFT. The Perdew-Zunger parametrization of local 
density approximation (LDA) (Perdew & Zunger 1981) is used for the exchange-correlation 
potential. Similar to the calculation of Hydrogen at the density of 80 g/cm 3 (Dai et al. 
2010; Recoules et al. 2009), the Coulombic pseudopotentials with a cutoff radius of 0.005 
a.u. for H and He are adopted. The plane wave cutoff energy is set to be from 200 Ry to 400 
Ry with the increase of temperature. As discussed in Ref. (Recoules et al. 2009), pressure 
derealization promises that the upper band electronic eigenstates are nearly plane waves, 
and thus the basis set is greatly reduced. Therefore, even though we choose a large number 
of bands in order to ensure the accuracy of the calculation, the computational cost is not 
very expensive. In our cases, we used enough bands in order to make the corresponding 
band energies higher than 8/c#T (especially at high temperature). Gamma point only is 
used for the representation of the Brillouin zone. We tested all these parameters carefully, 
and found that more k-points, more bands, larger energy cutoff, and more particles do 
not give any significant difference. The time step of QLMD is aj/ (20 y/kBTMj) (Recoules 
et al. 2009), where ai and Mj are the average ionic radius and the ionic mass of Ith ion, 
respectively. After thermalization, each structure is simulated for 5000 to 10000 time steps 
to pick up the useful information. 



3.2. Equation of state 

Firstly, we study the isochoric heating curve along the density of 10 g/cm 3 from 1 
eV to 550 eV under the conditions of the solar radiative zone. Comparisons between the 
pressure calculated using QLMD and AAB models for different compositions are given in 
Fig. 2. It is shown that for the system of He with two electrons at lower temperatures, the 
pressure of AAB model is much larger than that of QLMD. However, the relative difference 
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for H is not as big as for He. Therefore, with the increase of temperature and abundance of 
Hydrogen, the gap between QLMD and AAB models becomes less and less. It is very clear 
that more electrons and low temperature result in difficulties of describing the ionization 
balance for AAB model, where the composition of the system is complicated. The strong 
coupling among the ions, which is included in QLMD model but not in AAB model, plays 
an important role too in the calculations. Here, we recall that the simple way to calculate 
the EOS containing only one element and averaging the pressure according to the mass 
ratio is not accurate enough, as discussed in Ref. (Yuan 2002). For example, for the 
density-temperature point (10 g/cm 3 , 10 eV) and composition (X=0.7, Y=0.3, Z=0), the 
average pressure is P a = (0.7 x 301.3 + 0.3 x 58.6) = 228.49 Mbar. However, the pressure 
of this mixture obtained using QLMD is 214.77 Mbar, and the error is about 6%, which is 
significant for astrophysical applications. 

When the conditions go further into the solar interior, some interesting characteristics 
appear. For the density of 100 g/cm 3 , pressure-temperature relation is shown in Fig. 3. 
Moreover, the maximum density of the sun is about 160 g/cm 3 , whose isochoric heating 
curve from 10 eV to 1500 eV is shown in Fig. 4. Contrary to the relative difference between 
QLMD and AAB models in Fig. 2, the pressures of AAB model for 100 and 160 g/cm 3 are 
smaller than those of QLMD model, especially for one component H. In fact, the spherical 
assumption (Yuan 2002) for ions in AAB model makes the effective volume of the system 
larger than the real one, giving rise to smaller pressure. This phenomenon is obvious when 
the density is high enough, such as 160 g/cm 3 for H here. With the increase of temperature, 
the pressure of AAB and QLMD becomes consistent. At high temperatures, the ions are 
weakly coupled and electrons are almost free. Therefore, the semiclassical methods such as 
Thomas-Fermi and AAB can work well. Furthermore, it can be deduced that the EOS of 
the whole area of the sun can be obtained accurately and smoothly, which is beyond the 
OP and OPAL models, and completely parameter-free. 
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For the central conditions in the sun, the density reaches up to 152.7 g/cm 3 and 
the temperature up to 15696000 K according to the standard solar model (Bahcall et al. 
2001). The central pressure is also about 2.342xl0 5 Mbar, which is never reached from 
first principles simulation before. Based on QLMD, we solved this problem from an ab 
initio approach and obtained the EOS under this condition. We firstly compared the 
pressure close to the solar center with other models, as shown in Table. 1, and they agree 
reasonably well with each other. It is interesting to find that the pressure of AA model 
without energy- level broadening (AANB) is much larger than those of QLMD and AAB, 
and consistent with other models of EFF, Livermore (LIV) and MHD (Dappen et al. 1990). 
Since QLMD can be realized much more reliably, we can think AAB more reasonable, and 
both QLMD and AAB can improve the accuracy of EOS much. In fact, for the very dense 
medium, the electronic structures of elements are not only the energy levels, but also energy 
bands. It is reasonable to consider the effect of energy-level broadening, which is naturally 
included in QLMD and thus gives rise to appropriate results. For the solar cental regime, 
we calculate the pressure with different compositions. First of all, we consider the chemical 
compositions of H and He, with abundances X=0. 34828, Y=0. 65172, Z=0.0. It is found that 
the pressures of QLMD and AAB are almost equal within 0.5% error and are very similar 
to that of Standard. The small difference is caused by the lack of the heavier elements. In 
order to understand the effect of heavier elements on the EOS, we study the pressure of 
chemical compositions of H, He, C and H, He, O, respectively. As shown in Table. 1, the 
existence of heavier elements decreases the pressure a little, since the abundance is very 
small. In principle, we can calculate the conditions of solar interior with any composition 
and any abundance using QLMD resulting in EOS accurately and smoothly, particularly at 
relatively low temperatures and high densities. QLMD can also be very complementary to 
the AA models considering the computational cost, and improve significantly the accuracy 
of EOS in the solar interior. 
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In conclusion, the EOS of solar interior is calculated based on first principles method 
named QLMD, which is beyond the MHD and OPAL models. Furthermore, the model of 
QLMD can promise the accuracy for the EOS and be complementary to AA-like models. 
The accuracy can be improved significantly, particularly at relatively low temperatures and 
high densities where the coupling of ions can not be neglected. This gives us a useful tool 
to investigate the properties of the sun and sun-like stars from ab initio approaches, which 
is very powerful in astrophysics. This work can also open a new field to investigate the 
solar or stellar properties from first principles, promoting the study in astrophysics both in 
theory and experiment. 
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Table 1: Comparison of central pressures of the sun at typical temperature (T) and density 
(p). The results of EFF, LIV (OPAL) and MHD are in Ref. (Dappen et al. 1990) (X=0.34828, 
Y=0. 65172, Z=0). Results of AANB, AAB, and standard solar model (Standard) (Bahcall 
et al. 2001) are also compared. The pressures of chemical composition of H, He and C with 
X=0.3387, Y=6613, Z=0 ( a ), X=0.3125, Y=0.6406, Z=0.0469 ( 6 ), and chemical composition 
of H, He and O wi th X=0.3077, Y=0.6308, Z=0.0615 ( c ) are also shown. 

p (g/cm 3 ) T (eV) EOS model pressure (Mbar) 



141.25 


989.45 


EFF 


1.6396xl0 5 






LIV 


1.6162xl0 5 






MHD 


1.6190x10 s 






AANB 


1.6144x10 s 






AAB 


1.5922x10 s 






QLMD 


1.6090x10 s 


141.25 


1189.58 


EFF 


1.9579x10 s 






LIV 


1.9337x10 s 






MHD 


1.9365x10 s 






AANB 


1.9330x10 s 






AAB 


1.9027x10 s 






QLMD 


1.9122x10 s 


152.70 


1352.64 


Standard 


2.3420 xlO 5 






AANB a 


2.3603xl0 5 






AAB a 


2.3417xl0 5 






QLMD a 


2.3525xl0 5 






AANB 6 


2.2678xl0 5 






AAB 6 


2.2534xl0 5 






QLMD 6 


2.2575xl0 5 






AANB C 


2.2445xl0 5 






AAB C 


2.2301 xlO 5 






QLMD C 


2.2328xl0 5 
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Fig. 1. — (Color online) The coupling parameters (T) and degenerate parameters (6) of H 
and He from 1 eV to 550eV at 10 g/cm 3 (a)-(b), and from 10 eV to 1500 eV at the density 
of 160 g/cm 3 (c)-(d). 
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Fig. 2. — (Color online) Pressure vs temperature for the density of 10 g/cm 3 with different 
chemical compositions compared with results of AAB. In the figure, 0.7 represents X=0.7, 
0.4 represents X=0.4, respectively. 
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Fig. 3. — (Color online) Pressure vs temperature for the density of 100 g/cm 3 with different 
chemical compositions compared with results of AAB. 
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Fig. 4. — (Color online) Pressure vs temperature for the density of 160 g/cm 3 with different 
chemical compositions compared with results of AAB. 



